STABILITY OF VISCOUS SHOCKS IN ISENTROPIC 
GAS DYNAMICS 



BLAKE BARKER, JEFFREY HUMPHERYS 
KEITH RUDD, AND KEVIN ZUMBRUN 



Abstract. In this paper, we examine the stabihty problem for viscous 
shock solutions of the isentropic compressible Navier-Stokes equations, 
or p-system with real viscosity. We first revisit the work of Matsumura 
and Nishihara, extending the known parameter regime for which small- 
amplitude viscous shocks are provably spectrally stable by an optimized 
version of their original argument. Next, using a novel spectral energy 
estimate, we show that there are no purely real unstable eigenvalues for 
any shock strength. 

By related estimates, we show that unstable eigenvalues are con- 
fined to a bounded region independent of shock strength. Then through 
an extensive numerical Evans function study, we show that there is no 
unstable spectrum in the entire right-half plane, thus demonstrating nu- 
merically that large-amplitude shocks are spectrally stable up to Mach 
number M ~ 3000 for 1 < 7 < 3. This strongly suggests that shocks 
are stable independent of amplitude and the adiabatic constant 7. We 
complete our study by showing that finite-difference simulations of per- 
turbed large-amplitude shocks converge to a translate of the original 
shock wave, as expected. 



Consider the isentropic compressible Navier-Stokes equations in one spa- 
tial dimension expressed in Lagrangian coordinates, also known as the p- 
system with real viscosity: 



where, physically, v is the specific volume, u is the velocity, and p{v) is the 
pressure law. We assume that p{v) is adiabatic, that is, a 7-law gas satisfying 
p{v) = aoV~^ for some constants oq > and 7 > 1. In the thermodynamical 
rarified gas approximation, 7 > 1 is the average over constituent particles 
of 7 = {n + 2)/n, where n is the number of internal degrees of freedom of an 
individual particle [3]: n = 3 (7 = 1.66...) for monatomic, n = 5 (7 = 1.4) 
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for diatomic gas. More generally, 7 is usually taken within 1 < 7 < 3 in 
physical approximations of gas- or fluid-dynamical flow |27j . 

This system is an important and widely studied gas-dynamical model 
(see for example [27] and references within) , and yet little is presently known 
about the stability of its large-amplitude viscous shock wave solutions. Over 
two decades ago, Matsumura and Nishihara [24] used a clever energy esti- 
mate to show that small-amplitude shocks are stable under zero-mass pertur- 
bations. The linear portion of their analysis is equivalent to proving spectral 
stability, which through the more recent work of Zumbrun and collaborators 
[32l [2TI [23l [22I [3H [T7] , implies asymptotic orbital stability, hereafter called 
nonlinear stability. We remark that the results of (231 122| [31] hold for shocks 
of arbitrary amplitude, and thus nonlinear stability is implied by spectral 
stability. Hence, for large-amplitude shocks, spectral stability is the missing 
piece of the stability puzzle and the subject of our present focus. 

In this paper, we first improve upon the work in [24j slightly by extending 
the known parameter regime for which small-amplitude viscous shocks are 
provably spectrally stable, using an optimized version of the same method. 
We also show that this method cannot be extended any further to larger 
amplitudes. Using a novel spectral energy estimate, however, we are able 
to show that there are no purely real unstable eigenvalues for any shock 
strength. We note that this result is stronger than that which could be 
given by the Evans function stability index (sometimes called the orientation 
index), which only measures the parity of unstable real eigenvalues, see 
[HI HI ESI [3l|. A consequence of this result (which follows also by the 
Evans function computations used to determine the stability index [3T]) is 
that if an instability were to occur for large- amplitude viscous shocks, its 
onset, or indeed any change in the number of unstable eigenvalues, would be 
associated with a Hopf-like bifurcation in which one or more conjugate pairs 
of eigenvalues cross through the imaginary axis; see [281 [29| [30j for further 
discussion of this scenario. 

Continuing our investigations, we appeal to numerical Evans function 
computation to explore the large-amplitude regime through the use of wind- 
ing number calculations via the argument principle. Before doing so, how- 
ever, we rule out high-frequency instability through the combination of two 
spectral energy estimates, showing that unstable eigenvalues are confined to 
a bounded region independent of shock amplitude. This reduces the problem 
to investigation, feasible by numerics, of a compact set. Then by checking 
the low- frequency regime by repeating several Evans function computations, 
we determine whether or not a particular viscous shock is spectrally stable. 
As a final verification, we use a finite-difference method to simulate per- 
turbed large-amplitude viscous shocks, and check whether they converge, as 
expected, to a translate of the original profile. 

Conclusions and results of numerical investigations: We carry out 
our numerical experiments far into the hypersonic shock regime, exploring up 
through Mach number M « 3000 for 1 < 7 < 3. Particular attention is given 
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to the monatomic and diatomic cases, 7 = 5/3 and 7 = 7/4, respectively. In 
all cases, our results are consistent with spectral stability (hence also linear 
and nonlinear stability [23l [22l El]). This strongly suggests that viscous 
shock profiles in an isentropic are spectrally stable independent of both 
amplitude and 7. Our bounds on the size of unstable eigenvalues, which are 
independent of shock strength, may be viewed as a first step in establishing 
such a result analytically. 

Extensions and open questions. The present study is not a numerical 
proof. However, as discussed in [6], it could be converted to numerical proof 
by the implementation of interval arithmetic and a posteriori error estimates 
for numerical solution of ODE. This would be an interesting direction for 
future investigation. A crucial step in carrying out numerical proof by in- 
terval arithmetic is by analytical estimates special to the problem at hand 
to sufficiently reduce the computational domain to make the computations 
feasible in realistic time. This we have carried out by the surprisingly strong 
estimates of Section [5] and Appendices [BHCI 

A second very interesting direction would be to establish stability in the 
large-amplitude limit by a singular-perturbation analysis, an avenue we in- 
tend to follow in future work. Together, these two projects would give a 
complete, rigorous proof of stability for arbitrary-amplitude viscous shock 
solutions of ([1]) on the physical range 7 G [1,3]. 

2. Background 

By a viscous shock profile of ([T]), we mean a traveling wave solution 

v{x, t) = v{x — st), 
u{x, t) = u{x — st), 

moving with speed s and having asymptotically constant end-states {v^-, u±). 
As an alternative, we can translate x — > x — st, and consider instead sta- 
tionary solutions of 




Under the rescaling {x,t,v,u) — > {—esx, es'^t,v/e, —u/{£s)), where e is cho- 
sen so that < v-^ < V- = 1, our system takes the form 

Vt+Vx- = 0, 

ut + Ux + [av ^)x= { — , 

\ V /x 

where a = aoe~'^~^ s~'^ . Thus, the shock profiles of ([3]) satisfy the ordinary 
differential equation 

v' -u' = 0, 
u' + {av-"/)' = 




,,/\ ' 
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subject to the boundary conditions (v(±oo), M(iboo)) = {v±,u±). This sim- 
phfies to 

v' + {av-^y = (^^ 

By integrating from —cxd to x, we get the profile equation 
(4) v' = v{v -l + a{v-^ - I)), 

where a is found by setting x = +00, thus yielding the Rankine-Hugoniot 
condition 

Evidently, a 7"^ in the weak shock limit f+ 1, while a ~ in the 
strong shock limit w+ — > 0. 

Remark 2.1. Since the profile equation (|4]) is first order scalar, it has a 
monotone solution. Since f+ < V-, we have that Vx < for all x G M. 

By linearizing ([3]) about the profile {v,u), we get the eigenvalue problem 

Xv + v' — u = 0, 

(6) ^ , fh{v) 

V 



\u + u 

\V^^^ J \v 

where 

(7) h{v) = -v^+^ + a(7 - 1) + (a + l)v^ 

We say that a shock profile of ([1]) is spectrally stable if the linearized system 
([6]) has no spectrum in the closed deleted right half-plane 

P = {3f?e(A) > 0} \ {0}, 

that is, there are no growth or oscillatory modes for We remark that a 
traveling wave profile always has a zero-eigenvalue associated with its trans- 
lational invariance. This generally negates the possibility of good uniform 
bounds in energy estimates, and so we employ the standard technique (see 
[161 132j ) of transforming into integrated coordinates. This goes as follows: 

Suppose that (f , u) is an eigenfunction of ([6]) with a non-zero eigenvalue 
A. Then 

/X px 
u{z)dz, v{x) = / v{z)dz, 
-00 J —00 

and their derivatives decay exponentially as x ^ 00. Thus, by substituting 
and then integrating, {u,v) satisfies (suppressing the tilde) 

(8a) Xv + v' -u = 0, 

(8b) Xu + u' - ^:\{v' = — . 

This new eigenvalue problem differs spectrally from (l6|) only at A = 0, hence 
spectral stability of ([6]) is implied by spectral stability of ([8]). Moreover, 
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since ([8]) has no eigenvalue at A = 0, one can expect to have a better chance 
of developing a successful spectral energy method to prove stability. We 
demonstrate this in the following section. 



3. Stability of small-amplitude shocks 



In this section we further the work in [24] by extending slightly the known 
parameter regime for which small-amplitude viscous shocks are provably 
spectrally stable. We also show that this method cannot be extended any 
further for larger amplitudes. 

Theorem 3.1 ([2l])' Viscous shocks of ([T]) are spectrally stable whenever 



(9) 



,,7+1' 



07 



+ 2(7-1) 



07 / 



- (7 - 1) > 



In particular, as v^^- ^ 1 (hence 07 — > 1), the left-hand side of ([9|) ap- 
proaches 7 and so the inequality is satisfied. Therefore, small- amplitude 
viscous shocks of ([1]) are spectrally stable. 



Proof. We note that h{v) > 0. By multiplying (j8b|) by both the conjugate 
u and v'^'^^ /h{v) and integrating along x from cxd to —00, we have 



~h{d) 



+ 



u uv ' 



h{v) 



v'u 



h{v) 



Integrating the last three terms by parts and appropriately using 
substitute for u' in the third term gives us 



m to 



\\u\^v 



h{v) 



■ + 



h{v) 



+ / v{\v + v') + 



h{v) 



h{v) 



We take the real part and appropriately integrate by parts to get 



Ke(A) 



^7+1 



I |2 , I |2 

n + \v\ 



+ 



where 



g{v)\uY + 



+ 



h{v) 



\u 



l\2 



0, 



h(y) 



Thus, to prove stability it suffices to show that g{v) > on [^4., 1]. 
By straightforward computation, we obtain identities: 



(10) 
(11) 



'jh{v) — vh' {v) = 07(7 — 1) + i)'^^'^ and 
07 — h{v). 



u'u. 
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Using (jlOp and (llip . we abbreviate a few intermediate steps below: 



9{v) 



' 2 
' 2 



+ 



dv 



{-i + l)v^h{v)-v'^+^h'{i 

v"^ {{-f + l)h{v) - vh' {v)) d 
h(v)^ dv 



jv^ h{v) — v^h'{v) 
^h[v) — vh'{v) 



h{v) 



(a7 - h{v)) 



av^v 



7-1 



2h{v) 



[7^(7 + 1)^)^+' - 2(a + 1)7(7' - 1)*^+' + (a + 
+ 07(7 + 2)(72 - l)v - a{a + 1)7^(72 - 1)] 



(12) 



' 2h{vf 



[(7 + 1)1)^+2 + z)T(7 - 1) ((7 + l)v - (a + 1)7)^ 



> 



+ 07(7^ - 1)(7 + 2)v - a{a + 1)7^(7' - 1)] 

[(7 + l)i)^+2 + 07(72 - 1)(7 + 2)v - a{a + 1)7^(7' - 1)] 



(13) 



> 



2h{vf 

72a^t)^(7 + 1) 
2h{v)^v+ 



Thus from Q, we have spectral stability. 



□ 



□ 



We note that the hypothesis in ^ is not sharp. Indeed, one can show 
from (fT2]l . that a stronger condition could be given as 

(7 + l)t)^+2 + {)^(7 - 1) ((7 + l)v - (a + 1)7)2 



(14) 



+ 07(72 - 1)(7 + 2)v - a{a + 1)7^(7^ - 1) > 0, 



which is sharp in the following sense: When this condition fails to be true, 
then g{v) is no longer nonnegative, and thus the energy method fails. In 
Figure dja), we see that g{v) dips on the left-hand side when this inequality 
is compromised. We remark further that near vj^, the left-hand side of ([1] 
is monotone increasing in v. Thus, (jl4p holds if and only if 



(15) 



(7 + l)^'^' + ^^^(7 - 1) ((7 + 1)^^+ - (a + 1)7)' 
+ 07(72 - 1)(7 + 2)v+ - a[a + 1)72(72 - 1) > 0. 



We see from Figure [I{b) that ()15p is only a marginal improvement over 
([9]). However, since (fT5]l is sharp, we cannot hope to prove large-amplitude 
spectral stability using this approach. Instead, we proceed by a combined 
analytical and numerical approach as in (BJ [71 [8] , first showing that unstable 
eigenvalues can occur only in a bounded set, then searching for eigenvalues 
in this set by computing the Evans function numerically. Before doing so. 
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Figure 1. In (a), we have a graph of g{v) against v for 
v+ = 1 X 10^'^ and 7 = 2.0. Note that g{v) dips down below 
zero on the left-hand size. Hence, the energy estimate will 
not generalize beyond the small-amplitude regime. In (b) we 
graph the stability boundaries (dark lines) given by ([9]) and 
(jl5p . where 7 and v+ are the horizontal and vertical axes, 
respectively. We see that in our scaling (fT5]) is only a modest 
improvement over ([9]). The dotted lines correspond from top 
to bottom as the parameter regimes for Mach numbers 2, 5, 
and 10 (see the appendix to see how to determine the Mach 
number). Hence, this energy estimate does not even hold for 
shocks at Mach 2 and 7 > 1.084. 



however, we show in the following section that real unstable eigenvalues do 
not exist, even for large-amplitude viscous shocks. 

4. No REAL UNSTABLE EIGENVALUES 

In this section we use a novel spectral energy estimate to show that there 
are no purely real unstable eigenvalues for any shock strength. We note that 
this result is stronger than that which could be given by the Evans function 
stability index (sometimes called the orientation index) , which only measures 
the parity of unstable real eigenvalues, see [U [151 EI]- The fact that this 
holds for all shock strengths is interesting because it is among the strongest 
statements about large-amplitude spectral stability that has been proven to 
date. 

Theorem 4.1. Viscous shocks of ([T|) have no unstable real spectra. 

Proof. We multiply (pb|) by the conjugate v and integrate along x from 00 
to —00. This gives 

Jr Jk Jr v^'^ Jr V 
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Notice that on the real hne, A = A. Using (i8a|) to substitute for Xv in the 
first term and for u" in the last term, we get 



l{u' — v') + / u 



h{v)v'v f {Xv' + v")v 



Separating terms and simplifying gives 

, f , f h(v)v'v , f v'v f v"v 
uu' + 2 u'v- / = X /—+/—. 

We further simplify by substituting for u' in the second term and integrating 
the last terms by parts to give, 

, , , h(v)v'v , f v'v f Vx , '"'''^ 

uu' + 2 I {Xv + v')v - I =X I —+ I ^v'v 



v'^^'^ Jr V Jr V 



which yields 



uu' + 2Xl \v\' + 2 I v'v-a^ I ^+ !K = x! ^. 



v^^'^ Jr V Jr V 



By taking the real part, we arrive at 



l\2 



This is a contradiction when A > 0. □ □ 

We remark that the absence of positive real eigenvalues limits the admis- 
sible onset of instability to Hopf-like bifurcations where a pair of conjugate 
eigenvalues crosses the imaginary axis. In the following section, we give 
an upper bound on the spectral frequencies that are admissible. In other 
words, we show that if a pair of conjugate eigenvalues cross the imaginary 
axis, they must do so within these bounds. 



5. High-frequency bounds 

In this section, we prove high-frequency spectral bounds. This is an im- 
portant step because it gives a ceiling as to how far along both the imaginary 
and real axes that one must explore for spectrum when doing Evans function 
computations. Indeed if we want to check for roots of the Evans function in 
the unstable half-plane, say using the argument principle, then we need only 
compute within these bounds. If no roots are found therein, then we have 
strong numerical evidence that the given shock is spectrally stable. (Indeed, 
at the expense of further effort, such a calculation may be used as the basis 
of numerical proof, as described in [SI [7].) In this section we show that the 
high-frequency bounds are quite strong, only allowing unstable eigenvalues 
to persist in a relatively small triangle adjoining the origin. Moreover, these 
bounds are independent of the shock amplitude. 
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Lemma 5.1. The following identity holds for 5ReA > 0; 

1 



(5Re(A) + |9m(A)|) / v\u 



2 / + 



,'|2 



(16) 



< V2 



h{v)^. 



V \ \u\ + / v\u \u 



Proof. We multiply (|8bp by vu and integrate along x. This yields 



7*1 



u\ + vuu 



.'|2 



h(v) 

— V u. 



We get p6p by taking the real and imaginary parts and adding them to- 
gether, and noting that |Ke(z)| + 1 3=771(2;) | < \/2|2;|. □ □ 



Lemma 5.2. The following identity holds for ?R.eX > 0; 



(17) 



u 



n2 



23f?e(A)2 / |7;|2 + sRe(A) 



vf 1 

R V 2 



h{v) ^ a-f 



Proof. We multiply (|8bp by v' and integrate along x. This yields 



,.'|2 



A / uv' + / ti'u 



1 

R V 



u"v' 



1 



{Xv' + t;")?)'. 



Using ()8ap on the right-hand side, integrating by parts, and taking the real 
part gives 



3?e 



A / uv' + / uv 



h{v) 



+ 



{;7+l 2{)2 

The right hand side can be rewritten as 

1 



\v'\^ + ^e{X) [ 



m'|2 



(18) 5Re 



A / uv' + / u'v 



h{v) 07 



\v'\^ + ^e{X) I 
Jm 



,/|2 



Now we manipulate the left-hand side. Note that 

X [ uv' + [ u'v' = (A + A) [ uv' - [ u{Xv' + v") 



-23f?e(A) u'v- uu" 
Jr Jr 

-23f?e(A) / {Xv + v')v + 



.'1 2 



u 



Hence, by taking the real part we get 



A / uv' + / u'v 



This combines with (jlSp to give (jl7p . 



./|2 



25Re(A)^ 



□ 



□ 



Remark 5.3. Lemma \5.S\ is a special case of the high-frequency hounds 
given in [22l[31]. We also note that (fTTP follows from a "Kawashima-type" 
estimate as described in \22\ [3T] . 
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Lemma 5.4. For h{v) as in ([7]), we have 
(19) 



sup 

V 



1-V+ 

Ir. < 7- 

I — vl 



Proof. Defining 

(20) H(v) := h{v)v-"< = -v + 0(7 - l)v-"< + (a + 1), 

we have H'{v) = —1 — 07(7 — l)v < for < f+ < w < v_ = 1, hence 
the maximum of on t) E [ii-i-jW^] is achieved aX v = Vj^. Substituting (0) 
into (j20|) and simphfying yields (|19|) . □ 

We complete this section by proving our high-frequency bounds. 

Theorem 5.5. Any eigenvalue A 0/ ([5]) with nonnegative real part satisfies 

(21) 



1, 



3f?e(A) + |9m(A)|<(V7 + ^)'- 

Proof. Using Young's inequality twice on right-hand side of (jl6p together 
with (fT9|) . we get 



(5Re(A)+|9m(A)|) / vj-up 
Jr 



2 / ^a;!-"! + 



/|2 



< \/2 



< 



h{v) 
7)7+1 



(^2)2 /■ /i(t}) 



4^ 



-i}|ur + e / vlti'P + 



4e 



Assuming that < e < 1 and 1 



I '|2 I 

r" + 



7^1 
29 4e 



urn 



^1 — e)/2, this simplifies to 



(3f?e(A) + |9m(A)|) / + (1 - e) 

1 - e 



.'|2 



< 



7^1 
20 4e 



Applying (fT7|) yields 

(Ke(A) + |9m(A)|) / {)|n|2 < 
or equivalently, 

(3f?e(A) + |9m(A)|) < 
Setting e = 1/(2^+1) gives ([21]). 
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1 ■ 

1-e ^ 4^ 



(47 - l)e - 1 
4e(l-e) 



□ 



□ 



In the following section, we do an extensive Evans function study to ex- 
plore numerically the rest of the right-half plane in an effort to locate the 
presence of unstable complex eigenvalues. 
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6. Evans function computation 

In this section, we numerically compute the Evans function to locate 
any unstable eigenvalues, if they exist, in our system. The Evans function 
D{X) is analytic to the right of the essential spectrum and is defined as the 
Wronskian of decaying solutions of the eigenvalue equation for the linearized 
operator ([8]) (see [H [ini [HI ttSJ [I3l [l] ) . In a spirit similar to the characteristic 
polynomial, we have that -D(A) = if and only if A is in the point spectrum 
of the linearized operator ([8]). While the Evans function is generally too 
complex to compute explicitly, it can readily be computed numerically, even 
for large systems p9]. 

Since the Evans function is analytic in the region of interest, we can nu- 
merically compute the winding number in the right-half plane. This allows 
us to systematically locate roots (and hence eigenvalues) within. As a re- 
sult, spectral stability can be determined, and in the case of instability, one 
can produce bifurcation diagrams to illustrate and observe its onset. This 
approach was first used by Evans and Feroe [13] and has been advanced 
further in various directions (see for example [261 [ISl [H [SI [3 [5l [19] ) . 

We begin by writing ([8]) as a first-order system W' = A(x, X)W, where 



A 1 

(22) A{x,X) = I 1 



d 

dx ' 



Xv f{v)-X^ 

and f{v) = V — v^'^h{v), with h as in ([7|). Note that eigenvalues of ([8|) 
correspond to nontrivial solutions of W{x) for which the boundary condi- 
tions W{zizoo) = are satisfied. We remark that since v is asymptotically 
constant in x, then so is ^(x,A). Thus at each end-state, we have the 
constant-coefficient system 

(23) W' = A^{X)W. 

Hence solutions that satisfy the needed boundary condition must emerge 
from the unstable manifold Wi{x) at x = —oo and the stable manifold 
W^(x) A W^{x) at X = cxD. In other words, eigenvalues of correspond 
to the values of A for which these two manifolds intersect, or more precisely, 
when D{X) := det{W^W^W^)\^=o is zero. 

As an alternative, we consider the adjoint formulation of the Evans func- 
tion \26\ [1], where instead of computing the 2-dimensional stable manifold, 
we find the single trajectory W^{x) coming from the unstable manifold of 

(24) W' = -WA{x,X) 

at X = oo. Note that W^{x) is biorthogonal to both W2~{x) and W^{x) 
since {W{x)W{x)y = for all x and the initial data of W is biorthogo- 
nal to that of W. Hence, the original manifolds intersect when and 
are biorthogonal. Therefore, the adjoint Evans function takes the form 
D+{X) := (W^W,-)\.=o- 
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L 


7 = 1.2 


7 = 1.4 


7 = 1.666 


7 = 1.8 


8 


1.23(-1) 


1.16(-1) 


1.08(-1) 


1.04(-1) 


10 


2.07(-2) 


1.46(-2) 


1.75(-2) 


1.78(-2) 


12 


2.00(-3) 


1.40(-3) 


9.85(-4) 


7.20(-4) 


14 


6.90(-4) 


5.31(-4) 


4.73(-4) 


4.71(-4) 



Table 1. Relative errors in D{\) are computed by tak- 
ing the maximum relative error for 60 contour points eval- 
uated along the semicircle in Figure [2][^a). Samples were 
taken for varying L and 7, leaving vj^ fixed at = 10~^ 
(Mach M 1669). We used L = 8,10,12,14,16 and 
7 = 1.2, 1.4, 1.666, 1.8. Relative errors were computed us- 
ing the next value of L as the baseline. 



To further improve the numerical efficiency and accuracy of the shooting 
scheme, we rescale W and W to remove exponential growth/decay at infinity, 
and thus eliminate potential problems with stiffness. Specifically, we let 
VF(x) = ^V{x), where ^~ is the growth rate of the unstable manifold at 
X = —00, and we solve instead V'{x) = {A{x^ A) — fi~ I)V{x). We initialize 
V{x) at X = —00 to be the eigenvector of ^-(A) corresponding to 
Similarly, it is straightforward to rescale and initialize at x = 00. 

Numerically, we use a finite domain for x, replacing the end states x = 
±00 with X = zizL, for sufficiently large L. Some care needs to be taken, 
however, to make sure that we go out far enough to produce good results. 
In Appendices [B] and \C\ we explore the decay rates of the proffie v and 
A(x,\) and combine our analysis with numerical convergence experiments 
to conclude that our domain [—L, L] is sufficiently large. Our experiments, 
described below, were primarily conducted using L = 12, with relative error 
bounds ranging mostly between 10~^ and 10~^. Larger choices of L were 
needed on the high end of the Mach scale, going up to L = 18 in some cases, 
to get the relative errors down to 10~^. In Table [H we provide a sample of 
relative errors in D{X) for large-amplitude shocks. 

We remark also that in order to produce analytically varying Evans func- 
tion output, the initial data. V{—L) and V{L), must be chosen analytically. 
The method of Kato [20l pg. 99] , also described in [Sj , does this well by re- 
placing the eigenvectors of ([23j) with analytically defined spectral projectors 
(see also [SllIH]). 

Before we can compute the Evans function, we first need to compute 
the traveling wave profile. We use both Matlab's ode45 routine, which 
is the adaptive fourth-order Runge-Kutta-Fehlberg method (RKF45), and 
Matlab's bvp4c routine, which is an adaptive Lobatto quadrature scheme. 
Both methods work well and produce essentially equivalent profiles. 

Our experiments were carried out on the range 

(7,M) G [1,3] X [1.6,3000]. 
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(a) (b) 

Figure 2. The graph of the contour and its mapping via 
the Evans function. In (a), we have a contour, which is a 
large semi-circle aligned on the imaginary axis on the right- 
half plane (horizontal axis real, vertical axis imaginary). In 
(b) we have the image of the contour mapped by the Evans 
function. Note that the winding number of this graph is 
zero. Hence, there are no unstable eigenvalues in the semi- 
circle. Together with the high-frequency bounds, this implies 
spectral stability. Our computation was carried out for 7 = 
5/3 and u+ = 1 x 10~^. This corresponds to a monatomic 
gas with Mach number M ~ 1669 (see Appendix |A|) . 

Recall, for 7 E [1,3], that shocks are known to be stable for M S [1, 1.6], by 
([15]) , Section [3l hence this completes the study of the range 

(7,M) G [1,3] X [1,3000] 

from minimum Mach number M = 1 far into the hypersonic shock regime, 
and encompassing the entire physically relevant range of 7. 

For each value of 7, the Mach number M was varied on logarithmic scale 
with regular mesh from M = 1.6 up to M = 3000. We used 50 mesh points 
for 7 = 1.0 + O.l/c, where k = 1,2, .. . ,20. For the special values 7 = 1.4 
and 1.666 (monatomic and diatomic cases), we did a refined study with 400 
mesh points. 

For each value of (7,M), we computed the Evans function along semi- 
circular contours that contain the triangular region found in the previous 
section via our high-frequency bounds; see Figure El^a) . The ODE calcula- 
tions for individual A were carried out using Matlab's ode45 routine, which 
is the adaptive 4th-order Runge-Kutta-Fehlberg method (RKF45), and after 
scaling out the exponential decay rate of the constant-coefficient solution at 
spatial infinity, as described in [U [71 [U EJ [19] . This method is known to 
have excellent accuracy [H [19]; in addition, the adaptive refinement gives 
automatic error control. Typical runs involved roughly 300 mesh points, 
with error tolerance set to AbsTol = le-6 and RelTol = le-8. Values of 
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A were varied on a contour with 60 points. As a check on winding number 
accuracy, it was tested a posteriori that the change in argument of D for 
each step was less than 7r/25. Recall, by Rouche's Theorem, that accuracy 
is preserved so long as the argument varies by less than tt along each mesh 
interval. 

In all the cases that we examined, the winding number was zero. This 
indicates that the shocks we considered are spectrally stable, and in view 
of [2l| [23l I22j . nonlinear stability follows. Moreover, in light of the large 
Mach numbers considered, this is highly suggestive of stability for all shock 
strengths. 

7. Numerical evolution of strong shocks 

In this section, we use a standard finite-difference method to simulate 
perturbed large-amplitude viscous shocks, and show that they converge, as 
expected, to a translate of the original profile. 

We do this with a nonlinear Crank-Nicholson scheme with a Newton solver 
to deal with the nonlinearity. This method provides second-order accuracy 
and will allow for larger time steps than a naive explicit scheme. By differ- 
entiating the viscosity term in ([3]), we have 

Vt + Vx-Ux = 0, 



, --7— 1 "'XX 

Ut + Ux- ajv ' Vx = — 



V 

By implementing the Crank-Nicholson averaging, we obtain 
^-+1 - 1 



and 



^ _|_ __/'„"+l _ _|_ _ ^ 



At +4Ax^^J+i ""j-^' 



+ 16(Ax)2(^;p2("i^l - "i-l' + - ^i-l) 

,n+l „,n+l 



where n and j are, respectively, the discretized temporal and spacial indices. 
To cope with the nonlinearities, we use the Newton solver 



k 
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Figure 3. Snapshots of the evolution of a perturbed vis- 
cous shock wave solution generated by our extended Crank- 
Nicholson scheme (top curve: v against x; bottom curve: 
u against x, with time t fixed and increasing from fig- 
ure to figure). The parameters used were 7 = 1.4 and 
W4. = 9 X 10~^. This corresponds to a diatomic gas with 
Mach number M ^ 2877 (see Appendix |A]). As expected, 
the wave converges to a translate of the original shock. 



where F(C/"+\F"+1) and G(?7"+\F"+1) are the above finite difference 
schemes, and D^F, D^F, D^G, and D^G are their corresponding partial 
derivatives. Hence, we use the previous time step as our initial guess in 
Newton's method and then iterate until convergence. 

In Figure [3l we see the evolution of a perturbed viscous shock. We used 
a perturbation with a positive mass so that we could observe the conver- 
gence to a translate of the original profile, thus numerically demonstrating 
nonlinear stability. 



Appendix A. Mach number for the ^-system 
The Mach number is defined as 

u+ — a 
M = — , 
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where ti+ is the downwind velocity, a is the shock speed, and c+ is the down- 
wind shock speed, all in Eulerian coordinates. By considering the conserva- 
tion of mass equation, we have pt + {pu)x = 0. Hence, the jump condition 
is given by a[p\ = [pn], which implies, in the original scaling for ([T]), that 

U+V- — 

a = . 

V- — v+ 

Hence, 

^ ^ u+-a ^ v+{u- -u+) ^ v+[u] ^ _^v+ 
c+ c+{v--v+) c+[v] c+ 



or 



c+ J -p'{v+) 700^+'^ ^ 7 oo' 

To express this in our scaling, which is given in ([3]), we need to swap the 
pluses and minuses. Noting that < v+ < f_ = 1, we simplify to get 

M2 = — + = — . 

7f _j_ 1 — f + ja 
Recalling that a as — > 0, we find that 

(25) v+ ~ (7M2)-7 
as M — > 00. In particular, 

(26) |logt;+| ~7"^(log7 + 21ogAf). 

Appendix B. Profile bounds 
Denote profile equation ^ as v' = H{v,Vj^) := v{v — 1 -|- a{v~^ — 1)). 
Lemma B.l. For 7 > 1, < x < 1, 

27 1 < < 7. 

1 — X 

Proof. By convexity of f{x) = x'^, the difference quotient (j27p is increasing 
in X, bounded above by /'(I) = 7, and below by its value at x = 0. □ 

Proposition B.2. For 7 > 1 and < v < ^, < 

(28) --/{v-v+)<H{v,v+)<-^{v-v+). 
For 7 > 1, f+ < and j < v < V-. = 1, 

(29) ^{v-v^)<H{v,v+)<{v-v^). 
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Proof. By ([5]), 



H(v, va_) = v{ (v — I 



"((^--)+(}^)(er-0) 

{v-v+)(v 



Applying (j27p with x = we obtain (j28p from 

f — 7 < < V — (1 — va.) 

V — v+ 

Similarly, we may obtain ()29p by the calculation 



t> — 1 V \v/\l — v^/\l — V// 4 4 

□ 

Corollary B.3. For 7 > 1, < w+ < j^, and v{0) := f+ + the solution 
V of ^ satisfies 

(30a) -v+\ < (^^^e^^ x > 0, 

(30b) |{)(rE) - u_| < Q)e^ x < 0. 

Proof. Observing that {v — v^){0) = we obtain (ISOap by Proposition 
IB.2l and the comparison principle for first-order scalar ODE. Likewise, ()30bp 
follows by (j29p together with the observation that, by convexity of H, \H\ 
is bounded below by estimates obtained at v = v+ + and w = | of ^ and 
|, respectively, so that v traverses [v+ + |] over an x-interval of length 
< ^ - 12 □ 

- 1/16 ~ '-' 
Remark B.4. From Proposition \B.S\ and Corollary \B.^ we obtain the re- 
markable fact that, in the scaling we have chosen, v decays up to first de- 
rivative to endstates Vj- at a uniform exponential rate independent of shock 
strength, despite the apparent singularity at v = v+. 

Appendix C. Initialization error 

Lemma C.l. For A as in (j22p . | • | the Euclidean (i"^) matrix operator 
norm, 

(31a) |^(..A)-^+(A)|<(M±i±^!(l^i)!!±l),-¥, ,>o. 
(31b) M(..A)-^-(A)|<( ^l^l + '+^''^-" )e^. x<a 
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Proof. By ([22]), |>l(x,A)-^±(A)| < 2X\v-v±\ + \f{v)- f{v±)\. As computed 
in the proof of Lemma [5^ f'iv) = — 1 — 07(7 — l)v~^~^ . Applying (|27l) to 
the expression for a in ([5]), we obtain < a < 7^^, so that 

\f'{v)\ < 1 + 7^(7 - 1)^^"' < 1 + 7^(7 - 1)^^;', 

yielding (ISlah by the Mean Value Theorem and ()30ap . Bound (I31bp follows 
similarly. □ 

Theorem C.2 (Simplified Gap Lemma [151 [6l [7]). Let and Jr^ be a left 
eigenvector and associated eigenvalue of —A'^{\) and suppose that 

(32) ^ ' - 

|(A- A+)(x)| < Cze"''^ x>0, 

with < fj < r]. Then, there exists a solution W = e'^^^V{x, A) of ()24p with 

(33) \Vi^,X)-V^X)\ C,C.e-^^ 

\V+{X)\ iv-r]){l-e) 

provided (rj — fi)^^CiC2e~^^ < e. Similar estimates hold for solutions of 
(I22D onx<0. 



Proof. Writing V' = V{—A~^ — fi'^) + V{—A + A~^) and imposing the limiting 
behavior V{+oo,X) = , we obtain by Duhamel's Principle 



+00 



V{x) = V+- y(y)e(-^^-^^)(^-^)(^ - A+)dy, 



from which the result follows by a straightforward Contraction Mapping 
argument using ([32]) . See [IS El [7] for details. □ 

Lemma C.3. For A satisfying (I2ip . 7 G [1,3], and < j^, l\32\\ (a) holds 
with rj = 1/27, f) = 1/47, and Ci = 10^, and similarly for x < 0. 

Proof. Using the inverse Laplace transform representation 

g(-A+-M+)x = J_ / e^\z + + Jl+)-^dz, 
2m /r 

where L is a contour enclosing the eigenvalues of {—A"^ — pt+) and distance 
fj away, and estimating the resolvent norm \ {z + A'^ + Jl'^)^^\ by Kramer's 
rule, we obtain the stated crude bound, and similarly for x < 0. □ 

For error tolerance := (| log0| ~ 2k), define 

L_(0) := 2(1 log 10-^1 + I log(27 + 7 + 27=^(7 - 1))| + | log 0|) + 12, 

L+{e,v+):=^(\log 10-^1 + I log(27 + 7 + 7^(7 - l)z;;i)| + I log ( 
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Corollary C.4. For A satisfying (I2ip . we have relative error hounds 

|VFf(-L_,A)-y-e^-1 \Wt{L+A)-V^ei'*-\ ^ ^ 
|y-e/^~^| ' llZ+e^^+^'l 

Proof. This follows by (l33j) with ([31]), Lemma [Ol and ([21]). □ 

Remark C.5. Combining (j34p mi/i (|26p . we /inc? i/iai 

L+ ~ ^(2 log M + 4 + I log 10-^1 + I log 0|) 

suffices to obtain relative initialization error less than 9 for 7 G [1,3] and A 
in the computational region (j2ip . For = 10"'^, we thus obtain 9 tolerance 
up to M = 3, 000 ~ 10'^ for ~ 40. This is conservative, as we have made 
little effort to optimize bounds, hut still within the realm of our experiments. 
Our numerical convergence studies indicate that L± = 18 in fact suffices for 
10~^ accuracy. 
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